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Abstract 


In this article, we formulate an entropy stable scheme for nonlinear degenerate convection diffusion 
equation with viscous terms in non-conservative formulations. Here in we extend the S.Jerez idea of first 
order entropy stable convection scheme for semi discrete scheme in to fourth order scheme using finite 
central difference scheme. Major advantage of this work is fourth order accuracy of the solution and fixed 
numerical diffusion term which can provide the non- oscillatory solution. Finally, few computational 
analyses are given to shown the accuracy of entropy stable scheme for degenerate parabolic equations. 
Keywords: Entropy stable scheme, Hyperbolic conservation law, Convection diffusion equation, Central 
difference operator. 
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1. Introduction: 
In this section we discussed some basic definitions 
and results. 


We need derive fourth order accurate aproximation 
to first and second derivative with respect tox 














1.1 Notations of difference operator | [We Ow i he O = + O(h*), 
We are discussing some general notations and Ox 6 Ox (5) 
basic formulas from [9]. The central difference 2 a2 
h” 0°w. Ow 4 
with respect to ¢ are defined by Se 
6 Ox” Ox (6) 
wet - wt 2 
= — = (1) = 14 82) ™ + O(n’), 
2 6 Ox (7) 
Similarly, with respect tox is Where, 
Wit Wat 
= n 2 2D 
: 2 CM ee AOU), (8) 
Ox 
The second central difference is defined by 
O-=5'°6 . Otw 
_ 4 
Here we have yt pO), ) 
So 2 Win ~2W, + Wi 
anes ke (3) Rewrite above equation (7) by 
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Oia 


5°) '6,w+ Oh* 
Ax eo CO): 


using the Binomial expansion 


2 
Ow a 82)5, w+ 0(h"), (10) 
Ox 6 

Our main motivation is to introduced a fourth 
order non-oscillatory entropy stable scheme for 
non-degenerate convection diffusion equation. 
Next section 1.1 discussed the basic results of 
fourth order scheme. The present article 1s 
organized in the following way. Some preliminary 
results and definitions are given in sections 1.1 and 
2. The fourth order entropy stable scheme and new 
limiter are discussed in sections 3 and 3.1. 
Numerical results of proposed framework 
discussed in section 3.1 and 4. 

1.2 Forth order scheme for degenerate equation 
Let us consider the following equation (11) with 
non-conservative diffusion term in_- one 
Dimensional space. 


w, +f(w), = (k(w)w,),,(% DE RxR* 


W(x,0) = w,(X). (11) 
Where fbe a nonlinear flux function of conserved 
vector quantity 

w(x, t):RxR* >QeER" and k(w)EeR™ is a 
positive semidefinite diffusion matrix defined in 
Q. The diffusion term vanish in some sub spatial 
interval. For simplicity we denote w = w(x,f). 
This type of system represented by in porous 
media flow and two phase flow model [1-5]. 


If = k then equation (11) will be the following 


conservative form, 


w, +f(w), = K(w),,. (12) 


Due to this degeneracy of viscous term near shock 
may have solution in_ hyperbolic-parabolic 
dynamic problem. 

Let us discretized the equation (11) using central 
difference operator, 


Ne 
of (w) = pall Talla + O(h*) 
Ox Ox (13) 


Using Taylor series expansion, approximate first 
derivative by fourth order approximation in h. 


Similarly, 
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0(k(1+ — 657 )d,(w)(x) 
Ok(w)w, — 6 + O(h*), (14) 
Ox Ox 
By the product rule, R.H.S of the (14) is the 
following form. Let us consider 
ok(w)w,. _ HkL) +0(hA) 
Ox Ox 


Aen = (KL) (KL) 3.) + 


2h 12 pe 
4 > 4 
ae jt —(KL) j-l ns +O}h ) 
= I 1 1 
= 5 kad kj1)+O(h* ) C5) 
where 
1 -1 7 
K jas = A 1a wise Wie Lr oo Wat —Ww,) 
7 1 
oe (wW ~Wi2) Kj Wj-2 —W;3) 


7 
» (+k allwll . 


= k...[[w 
he i2 jee pall M3 


1 
+S kylbel gg blll, --> (16) 
We replaced the value of L by first order forward 
difference operator for representing the first 
derivative. Substitute the equation (16) in to (14), 
then we get a fourth order scheme. Replace the 
notation 


h = Ax. The numerical aa will be 





O Of (w(x, t Ow(x,t 
“ee - oa Ds Skok t) a7) 
Ot 
ea | 4 4 1 
“a Pave I Ax A (k jet K j- 
= 7 
Fjasy2 = (“qahine tha th 


~ 75 ita 
2. Entropy stable scheme for degenerate 
convection diffusion equation 
2.1.Basic results 
Let us consider (12) with K — 0, then the equation 
have hyperbolic in nature. The main challenge for 
solving such type of equation is the shock wave 
formation near discontinuities of solution profile 


al 
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of nonlinear flux. Even for smooth initial data, the 
solution profile may have shock formation with 
nonlinear flux function. So we can use the weak 
formation in that situation. It is well known that 
weak solution do not converge to physically 
relevant unique solution. Lack of uniqueness of 
numerical solution must be enforced additional 
criteria, we call it as entropy condition. Next we 
discuss the entropy stability of the scheme (18). 
The accuracy and efficiency of entropy stable 
methods and higher order of accuracy discussed in 
(2, 3, 4] and discontinuous Galerkin methods in 
[5S]. The entropy inequality for a degenerate 
parabolic problem Our main motivation of the 
study is to introduced the higher order entropy 
stable scheme. Entropy condition for system of 
equations (12) is defined in [6] Consider the three 
tuple (n,q,r) of functions from the set Q to R, with 
strictly convex such that 


dw = Nwhw> (19) 
dk 

hy = Nw os (20) 

WT ww — = 0,w € 2, (21) 


where the subscript w denotes the Hessian 
function. Recall that 7 = y(w) and g = q(w) are the 
entropy flux function, respectively, and the new 
function r = r(w) 1s named the diffusion entropy 
flux. Since 7 is strictly convex and v = nyis the 
entropy variable. Silvia Jerez Et.al. prove that 
degenerate form equation (12) is satisfy 


NW) + QW)x— FW) xx S 9, (22) 


The equation(12) said to be entropy stable if there 
exist two numerical entropy fluxes Qj: and’ i+1/2 : 
consistent with g and r,which satisfies the 
following condition (see [6]) 


d 1 

Te)O +7 (92-84) 
1 ! 
a(t anna) =v ey 
2.1.1. Entropy stable scheme with Entropy 
conservative flux 


Let us consider the Lefloch formula for higher 
order flux, 


2 
Fru = ye se oe 6 F (Wj—s,Wj—s+r) (24). 
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Theorem 2.1. (seel LSI, Theorem 4.4]). For p €N, 


assume that ++ solve the p linear equations 
p p 
» ra? 7 1» (25-1 (s=2,..9, ) (25) 
r=1 i=1 


then F‘2P}th-order accurate, in the sense that for 
sufficiently smooth solution w and_ entropy 
conservative by following discrete entropy 
equality. 

Oia y2 = Dh a DE=0 OCs Mj -str), (26) 
Substitute p = 2 in entropy conservative 

flux 


Z r—1 
ah _ p _— 
Prt /2 -_ » a, (Wj-s,Wj—s+r) 
r=1 s=0 
(27) 
and 
= f(w) + fs) 
Which imply that (28) 
1 1 
F jaaya = (-SRi2t — Fat os —F; asim tes 
(28) 
The two-point finite difference scheme, 
O lL : ~ 
ae = ee 
(29) 


~ 


Re (K lw) jatr2 - 


ky flwH 12) 


is entropy stable provided the numerical viscosity 
matrixk; 41/2 satisfies the following conditions 


Bia pakjasye Wi 4 = Ir]j412- (0) 


And[v]; 44/24 +1/2[W] 41/2 = 0. G1) 


proved in [6]. Here Jerez used _ entropy 
conservative flux, but due to the diffusion term 
numerical scheme become entropy stable provided 
k = 0.Suppose if k = O then scheme (29) will not be 
entropy stable. It satisfies only entropy 
conservative condition. In that case _ near 
discontinuity solution profile may exhibit Gibb’s 
phenomenon near discontinuity region. To prevent 
this phenomenon, we need to add extra numerical 
diffusion term. Jerez modified the scheme for 
entropy stability that we discussed in next 
subsection. 


a2 
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2.1.2. Entropy stable scheme with Entropy 
stable flux 
Consider the fourth order entropy stable flux 
defined in [2]. If K = O the numerical scheme 
equation will be 
a 1 
5 W(x, t) = —— Fat — Fy_4/2),32) 
replace Fj: by entropy stable fluxF;4, 4 /2> 

~ 1 
Flap = hapa 5 ALlvl] 263) 
Here D is positive definite matrix. The scheme 
(32) using flux (33) is entropy stable by Cheng 
theorem 3.3 in [7] and results in [2]. 
Under the hypothesis of the theorem 3.2 in [6] the 


finite the entropy stability of the difference 
scheme, 


0 1 
wnt) =—Z(F41—F,2) 4 


1 
a (lw I], ,2 — [Iw] } ay , (34) 
Wheree>0 is proved by Jerez. 


3. Fourth order entropy stable scheme for non- 
degenerate convection diffusion equation 
We consider fourth order entropy stable scheme 
for non-conservative formulation of diffusion term. 
For generating entropy stable fourth order scheme 
we need a central forth order scheme for 
convection diffusion equation with non- 
conservative diffusion term. 
Substitute the flux (33) to the two point scheme 
(29) is modified to, 
0 1 4 4 
mnt) = 3. (B— Bis) 4 
T 


Vv; ~ ~ 

a Ue 4 — k, Je?) 
Remark 3.1. The entropy stable scheme (34) by 
Jerez we discussed, but in that scheme diffusion 
term not fixed for PDE, it may vary due to the 
varying . In this article we used sign stable fourth 
reconstruction for diffusion term. ie., sign({[v]]) = 
sign(<~< v >>) where 


1if0<x<l 
—1, other wise 


sign(x) = (36) 


Theorem 3.1. The scheme (35) is entropy stable if 
it satisfies following condition. The numerical flux 
F'412 consistent with f and satisfies 
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[ViainPivi2 = [@lisi12 and the numerical viscosity 
matrix kj verifies 


ree = [r']j41/2, G7) 


and 
[v]) 1k 412 2 0G8) 
2 


Proof, Multiply v; by both side of equation (12) 


Tiff | _F _oF ltR 
Yj = (F,3 F1)= Pa pe 
1 


——!T" a 
(ll, -2F, 2) (40) 
The numerical fluxQ; + 21s defined by 


ote eee reg ree ey 


_— ma —1 
Similarlyv;' kj 41/2 = ve Ors k, allw]], = 


o =| 7 


Lille 2 Shallvl) yp 


After adding extra numerical diffusion term the 
equation (35) will be 
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0 l 14 14 
ay wat) = aan os Ca _ Ft) 
1 /- : 
Th (Fe41 7 k;_1) (44) 
4 — fF = ; Vil. 
Where the flux Huet — “j+s gD i42 it Ni+d 1S 


entropy stable flux 


Theorem 3.2. A scheme with numerical flux 
equation (44) is entropy stable with numerical 
entropy flux given by 


ao Gi 41/2 + V4 41/2D) 41/2 << V >> 4172 (45) 
is entropy stable with entropy flux is defined by 


2 = Gap 7, 42D; 110), 4) (40 


an — Vi41/2kj41/2 (47) 
with K v oe B2llvyj+12. ; 
Existing tools for making higher order 
reconstruction 1s WENO-Z and JS-WENO e.t.c. 
But generally, it is not sign stable. Here we 
reconstruct. 

[[v]]j+12 by fifth order WENO-Z method[8], but it 
is not sign stable. So we introduced limiter for 
preserving the sign stable. 

3.1. Limiter for sign stability 


In this article for reconstructing 5-4 we used 
WENO-Z scheme, but it is not sign stable. So, we 
used by following way 
1,if0<sG)<1 
= 48 
OW) 0, other wise oo 
where s(j) 1s the numerical jump at x(J) ie., 


si) = + (49) 


Vinge 
provided v;,, — uF 0 otherwise s(z) = 0. 
Here we consider the general Courant-Friedrich- 
Lewy(CFL) condition for parabolic equations is 
CFL = max|f,| = + 2 max|k| =. <1, (50) 
4.Test problems 
Consider the equation 


wet (5) = Kw)exs (x0) € [-2,2] X [0,7] 


(56) 

where w € [0,00] and diffusion matrix is defined as 
Kw) = uw, and K° (w) = k(w). Based on the 
known entropy flux pair (7,2) for the 

burger ao Entropy 3-tuple (7(w), 2g(w),r(w)) 


= (4, 4, 3uw*) is satisfied equation 1.11 in [6], 
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where uw = 0.01. Let us consider entropy stable 
scheme for that an entropy conservative numerical 
flux is consider. Fi+1/2 = 6(Mj41 + wie Tw/) 
and a numerical viscosity matrix satisfying 
theorem 3.4 in [6] by 
Kj 41/2 
1, ey Wj +1=0 


= 4w? + Wi Wj41 + Ww 
; other wise 
3 W; + Wy44 


Initial condition is x €[-2,2] 


if—-O0O5<x<0.5 


wo(x) =} 
a 0, other wise 


(51) 


For the computation we considering N = 1200for 
reference solution also. Mainly the following 
numerical methods are used for comparison 

ES1 flux is entropy conservative, nonconservative 
discretization of source term, numerical scheme is 
29. ES2 flux is entropy stable, non-conservative 
discretization of source term using the numerical 
scheme is 34. ES3 Fourth order entropy stable 
scheme. All test problem using ES3 are SSP-RK4 
method used for solving differential equations. 
Numerical simulations of ESI and ES2 are 
obtained combining an_ entropystable — spatial 
discretization with a TVD-RK2 time stepping. The 
errors and convergence rate of ESI and ES2 are 
shown in [6-9]. It is first order accurate. If u — 0 
the PDE have hyperbolic in nature. So spurious 
oscillation in solution profile will produce the 
numerical scheme near discontinuity. From figure 
1 we can understand after adding extra diffusion, 
oscillation removed. Next we are discussing the 
numerical result for non-oscillatory higher order 
entropy stable scheme. 

5. Test problems for fourth order non- 
oscillatory entropy stable scheme 

Recalling the equation (12). If uw tends to zero in 
(58) tends to hyperbolic case. Entropy 
conservative methods capture appearance of 
propagation of shock wave correctly but may 
produce strong oscillations in near shock region. In 
order to reduce the oscillation in [6] add some 
viscosity term. This work useful to capture 
correctly the non-classical shocks due to parabolic 
and hyperbolic interaction. For detailed reference 
see [6]. 
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wie 
w, + (=) = K(w),xx, (x,t) € [—2,2] x [0, T] 


(52) 
Test 1 





initial 
ES1 
reference solution 








Figure 1: Solution for first order entropy conservative 


flux with 100 points with initial condition(51) 
Example 5.1. [6] Let us consider the equation 
wA(w), = (K(w)w2)»(%t) € [-1,1]x[0, 1],w € Q 


(33) 
where Q = [0,+co]. 
with initial condition is defined by, 
(1 —x*)2, if-1<x<1 
= 4 
wolX) 0, other wise Ca 


,X € [—2,2]. Entropy conservative numerical flux 
is defined by 


, Solution at T=0.5 


~ — initial condiction 


© EC3 
0.8 


———— far 
0.6 


0.4 


0.2 


-1 1.5 


=e ee ee 
Me 
-1 5 
. 


-0.5 0.5 


condition (54) and N=100, u=0.01, CFL=0.9 
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where w € [0,00] and diffusion matrix is defined 


as K(w) = uw’, and K°(w) = k(w). 


Test 2 


1.2 
initial 
O ES2 

1 ———<—< reference solution 


0.8 


> 0.6 


0.4 


0.2 





45 2 


Figure 2: Solution for first order entropy stable flux 
with 100 points with initial condition(5 1) 


4 _F ,_1 
ie aaa Dally. (95) 
where D; +3 is maximum of first derivative of f and 


Vv is entropy variable. Numerical simulations of 
ES3 are obtained combining an _ entropystable 
spatial discretization with a TVD-RK4 time 
stepping. For that, we compute a reference solution 
with uw = 0.1 using the ES3 scheme with N = 1200 
grid values. 


Solution at T=1 


~~ initial condiction 
O EC3 
— ref 





x xX 
Figure 3: Numerical solution using ES3 with initial Figure 4: Numerical solution using ES3 with 


initial condition (60) and N=100, u=0.01, 
CFL=0.9 
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Example 5.2. Here we consider the equation 


Ww 


wet (5) = K(w)ex, (wt) € [-2,2] x 
[0, T | (56) 


with initial condition 














initial condiction 
EC3 
—™ 





0.8 


0.6 


0.4 


0.2 
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Wo (x) 
= if -O5<x<05 
= 0, 


other wise 


~ initial condiction 
© EC3 
— et 


Figure 5: Numerical solution using ES3 with Figure 6: Numerical solution using ES3 with 


initial condition (54) and N=200, u=0.01, CFL=0.5 


We used ES3 scheme with CFL=0.5, and Final 

time T= 1; Rate of convergence and accuracy table 

of the problem 57 is given. 

Conclusion 

Herein we introduce the concept of fourth order 

entropy stable scheme for degenerate convection 

equation and new limiter for preserving sign 

stability for numerical diffusion term. Using this 

limiter we can use all type of higher order WENO 

for reconstruction of diffusion term, Because of 

this we can generate a higher order nonoscillatory 

entropy stable scheme. 
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